{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from typing import Union\n",
    "\n",
    "import numpy as np\n",
    "from itertools import *\n",
    "from functools import *\n",
    "import operator\n",
    "import joblib\n",
    "import datetime\n",
    "from scipy.stats import norm\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"error\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "outputs": [],
   "source": [
    "# Signal Function\n",
    "def I(x: np.ndarray):\n",
    "    t = np.zeros([x.shape[0], 2])\n",
    "    t[:, 0] = x\n",
    "    return np.max(t, axis=1)\n",
    "\n",
    "# Expectation Function\n",
    "def E(x: Union[np.ndarray, list], p: np.ndarray):\n",
    "    return np.dot(x, p)\n",
    "\n",
    "# PDF function\n",
    "def p_norm(x: np.ndarray, mu, sigma):\n",
    "    init_arr = norm.cdf(x, mu, sigma)\n",
    "    res_arr = np.zeros(x.shape[0])\n",
    "    res_arr[0] = init_arr[0]\n",
    "    # res_arr[0] = init_arr[0] - norm.cdf(x[0] - 1, mu, sigma)\n",
    "    res_arr[1:] = np.diff(init_arr)\n",
    "    res_arr[-1] = 1 - init_arr[-2]\n",
    "    return res_arr\n",
    "\n",
    "# Variable\n",
    "alpha = 0.95\n",
    "# Kt\n",
    "mu_K_t = 6  # 均值\n",
    "cv_K_t = 0.3  # 变异系数\n",
    "sigma_K_t = cv_K_t * mu_K_t  # 标准差\n",
    "K_t = np.array([4, 5, 6, 7, 8])\n",
    "p_K_t = p_norm(K_t, mu_K_t, sigma_K_t)\n",
    "# Dt\n",
    "mu_t = 5\n",
    "sigma_t = 0.5\n",
    "D_t = np.array(list(range(0, 10)))\n",
    "p_D_t = p_norm(D_t, mu_t, sigma_t)\n",
    "# Joint Distribution of Kt, Dt\n",
    "K_D_t = [(K_t_i, D_t_i) for K_t_i, D_t_i in product(K_t, D_t)]\n",
    "# Consider as independent\n",
    "p_K_D_t = [p_K_t_i * p_D_t_i for p_K_t_i, p_D_t_i in product(p_K_t, p_D_t)]\n",
    "\n",
    "p = 19\n",
    "b = 19\n",
    "c_e = 6\n",
    "c_m = 2\n",
    "h_e = 1\n",
    "h_m = 0.2\n",
    "theta = h_m / (1 - alpha)\n",
    "\n",
    "max_epoch = 2\n",
    "z = range(-(max_epoch + 1) * np.max(D_t), (max_epoch + 1) * np.max(K_t) + 1)\n",
    "y = range(-(max_epoch + 1) * np.max(D_t), (max_epoch + 1) * np.max(K_t) + 1)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "outputs": [],
   "source": [
    "def g_1t(x, epoch):\n",
    "    val = h_e * E(I(x - D_t), p_D_t) + b * E(I(D_t - x), p_D_t) + (c_e - h_m) * x\n",
    "    if epoch != max_epoch:\n",
    "        # last cache need to be added for f_0t\n",
    "        res = alpha * E([f_0t(x - D_t_i, epoch + 1) for D_t_i in D_t], p_D_t)\n",
    "        val += res\n",
    "    return val\n",
    "\n",
    "def g_2t(x, epoch):\n",
    "    val = 0\n",
    "    if epoch != max_epoch:\n",
    "        # last cache need to be added for f_1t\n",
    "        val += alpha * E([f_1t(x - D_t_i, epoch + 1) for D_t_i in D_t], p_D_t)\n",
    "    return val\n",
    "\n",
    "def g_3t(x, epoch):\n",
    "    exp = [np.min([np.max([S_t(epoch), x]), np.max([x + K_t_i])]) for K_t_i in K_t]\n",
    "    val = (h_m - c_m) * x + c_m * E(exp, p_K_t)\n",
    "    if epoch != max_epoch:\n",
    "        val += (alpha * E([f_2t(x - D_t_i, epoch + 1) for D_t_i in D_t], p_D_t) +\n",
    "                alpha * E([f_3t(np.min([np.max([S_t(epoch), x]), x + K_t_i]) - D_t_i, epoch + 1) for K_t_i, D_t_i in K_D_t], p_K_D_t))\n",
    "    return val\n",
    "\n",
    "\n",
    "def f_0t(v_0, epoch):\n",
    "    return - c_e * v_0\n",
    "\n",
    "def f_1t(v_1, epoch):\n",
    "    return g_1t(np.max([v_1, S_2t(epoch)]), epoch)\n",
    "\n",
    "def f_2t(v_2, epoch):\n",
    "    return g_1t(np.max([v_2, S_2t(epoch)]), epoch) + g_2t(np.max([v_2, S_3t(epoch)]), epoch) - g_1t(np.max([v_2, S_2t(epoch)]), epoch)\n",
    "\n",
    "def f_3t(v_3, epoch):\n",
    "    return g_1t(np.max([v_3, S_4t(epoch)]), epoch) + g_2t(np.max([v_3, S_4t(epoch)]), epoch) + g_3t(np.max([v_3, S_4t(epoch)]), epoch) \\\n",
    "           - g_1t(np.max([v_3, S_3t(epoch)]), epoch) - g_2t(np.max([v_3, S_3t(epoch)]), epoch)\n",
    "\n",
    "\n",
    "# S_t can be cached\n",
    "S_t_cache = {}\n",
    "S_1t_cache = {}\n",
    "S_2t_cache = {}\n",
    "S_3t_cache = {}\n",
    "S_4t_cache = {}\n",
    "\n",
    "# init cache\n",
    "for c in range(1, max_epoch + 2):\n",
    "    S_t_cache.update({c : {\"value\": 0, \"cached\": False}})\n",
    "    S_1t_cache.update({c: {\"value\": 100000000, \"cached\": True}})\n",
    "    S_2t_cache.update({c : {\"value\": 0, \"cached\": False}})\n",
    "    S_3t_cache.update({c : {\"value\": 0, \"cached\": False}})\n",
    "    S_4t_cache.update({c : {\"value\": -100000000, \"cached\": True}})\n",
    "\n",
    "def S_t(epoch):\n",
    "    if not S_t_cache[epoch][\"cached\"]:\n",
    "        if epoch == max_epoch:\n",
    "            val = z[np.argmin([c_m * _z] for _z in z)]\n",
    "        else:\n",
    "            val = z[np.argmin([c_m * _z + alpha * E([f_3t(_z - D_t_i, epoch + 1) for D_t_i in D_t], p_D_t) for _z in z])]\n",
    "        S_t_cache[epoch][\"cached\"] = True\n",
    "        S_t_cache[epoch][\"value\"] = val\n",
    "        return val\n",
    "    else:\n",
    "        return S_t_cache[epoch][\"value\"]\n",
    "\n",
    "def S_1t(epoch):\n",
    "    return S_1t_cache[epoch][\"value\"]\n",
    "\n",
    "\n",
    "def S_2t(epoch):\n",
    "    if not S_2t_cache[epoch][\"cached\"]:\n",
    "        val = y[np.argmin([g_1t(_y, epoch) for _y in y])]\n",
    "        S_2t_cache[epoch][\"cached\"] = True\n",
    "        S_2t_cache[epoch][\"value\"] = val\n",
    "        return val\n",
    "    else:\n",
    "        return S_2t_cache[epoch][\"value\"]\n",
    "\n",
    "\n",
    "def S_3t(epoch):\n",
    "    if not S_3t_cache[epoch][\"cached\"]:\n",
    "        val = y[np.argmin([g_1t(_y, epoch) + g_2t(_y, epoch) for _y in y])]\n",
    "        S_3t_cache[epoch][\"cached\"] = True\n",
    "        S_3t_cache[epoch][\"value\"] = val\n",
    "        return val\n",
    "    else:\n",
    "        return S_3t_cache[epoch][\"value\"]\n",
    "\n",
    "\n",
    "def S_4t(epoch):\n",
    "    return S_4t_cache[epoch][\"value\"]\n",
    "\n",
    "\n",
    "def f_t(v_0, v_1, v_2, v_3, epoch):\n",
    "    return f_0t(v_0, epoch) + f_1t(v_1, epoch) + f_2t(v_2, epoch) + f_3t(v_3, epoch)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "outputs": [],
   "source": [
    "epoch = 1\n",
    "v_0 = 1\n",
    "v_1 = 4\n",
    "v_2 = 6\n",
    "v_3 = 8"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "outputs": [
    {
     "data": {
      "text/plain": "71.31434036663416"
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_t(v_0, v_1, v_2, v_3, epoch)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "outputs": [],
   "source": [
    "def C(t):\n",
    "    return c_e + np.sum([np.power(alpha, max_epoch + 1 - j) * h_e for j in range(t, max_epoch + 2)])"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C(1): 8.8525, C(2): 7.95, C(3): 7.0\n"
     ]
    }
   ],
   "source": [
    "print(f\"C(1): {C(1)}, C(2): {C(2)}, C(3): {C(3)}\")"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "outputs": [],
   "source": [
    "def relative_error(m, n):\n",
    "    return (m - n) / n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "outputs": [
    {
     "data": {
      "text/plain": "3.999999999999997"
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h_m / (1 - alpha)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}